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The reported new algorithm determines the exact exchange potential v x in a iterative way using 
energy and orbital shifts (ES, OS) obtained — with finite-difference formulas — from the solutions 
(occupied orbitals and their energies) of the Hartree-Fock-like equation and the Kohn-Sham-likc 
equation, the former used for the initial approximation to v x and the latter — for increments 
of ES and OS due to subsequent changes of v x . Thus, solution of the differential equations for 
OS, used by Kummel and Perdew (KP) [Phys. Rev. Lett. 90, 043004 (2003)], is avoided. The 
iterated exchange potential, expressed in terms of ES and OS, is improved by modifying ES at 
odd iteration steps and OS at even steps. The modification formulas are related to the OEP 
equation (satisfied at convergence) written as the condition of vanishing density shift (DS) — they 
are obtained, respectively, by enforcing its satisfaction through corrections to approximate OS and by 
determining optimal ES that minimize the DS norm. The proposed method, successfully tested for 
several closed-(sub)shell atoms, from Be to Kr, within the DFT exchange-only approximation, proves 
highly efficient. The calculations using pseudospectral method for representing orbitals give iterative 
sequences of approximate exchange potentials (starting with the Krieger-Li-Iafrate approximation) 
that rapidly approach the exact « x so that, for Ne, Ar and Zn, the corresponding DS norm becomes 
less than 10 -6 after 13, 13 and 9 iteration steps for a given electron density. In self-consistent 
density calculations, orbital energies of 10 -4 Hartree accuracy are obtained for these atoms after, 
respectively, 9, 12 and 12 density iteration steps, each involving just 2 steps of v x iteration, while 
the accuracy limit of 10" 6 -10" 7 Hart ree is reached after 20 density iterations. 

PACS numbers: 31.15.E-, 31.15.eg, 31.15.cj 



I. INTRODUCTION 



Modern ab initio calculations for molecular and con- 
densed matter systems are most frequently based on the 
density functional theory (DFT), which gives the pre- 
scription for determining properties of an interacting TV- 
electron system (in its ground state) with the aid of the 
Kohn-Sham (KS) scheme describing a corresponding sys- 
tem of non-interacting particles ■ While the DFT 
approach simplifies enormously the description of many- 
electron systems, application of the KS scheme requires 
the knowledge of the exchange-correlation (xc) potential 
which is the part of the KS Hamiltonian that represents 
all complex correlations between electrons in the physical 
system (due to their fermionic character and Coulomb in- 
teractions) . The xc potential is defined as the functional 
derivative v xca (r) = 6E xc [n^, n\\/ 5n a {r) of the xc energy 
E xc with respect to the density n a (r) of electrons with 
spin a (=t; i) but its exact form remains unknown, so 
it is usually found using the local-density or generalized- 
gradient approximations (LDA, GGA) to the functional 
E X c[n-\-,ni\. Although this approach has proved to be 
highly successful in calculations of numerous physical 
properties, it fails in some cases (e.g., for energies of 
bound unoccupied states) due to the incomplete cance- 
lation of Coulomb self-interaction in the KS potential 
including v xc<7 obtained with LDA or GGA. Such short- 
comings of these approximate xc potentials are largely re- 
moved when the exact exchange potential v xa (r), which 
is the dominant part of v XC(J = v xa + v ca , is used. It is de- 



fined as Vxo-(t) = SE x [na]/ '8n a (r) with the a-subsystem 
exchange energy E x given explicitly in terms of the oc- 
cupied KS orbitals which are obtained with the potential 
v sa that corresponds to n a (by virtue of the Hohenberg- 
Kohn theorem The KS potential v sa {r) including 

the so defined exact v xa is free from self-interaction and, 
consequently, it has the proper asymptotic dependence 
(— 1/r for atoms), unlike v sa (r) obtained with LDA or 
GGA. 

The KS orbitals <paa(f) and their energies t aa are 
eigensolutions of the KS equation 

hsa(r) (/)aa(r) = e aCT Q(T (r), C aa < € a +l,a , (1.1a) 

where 

h S a(r)=i(r)+v sa (r), t{r) = -\V 2 {r) (1.1b) 

(atomic units used throughout the paper, orbitals chosen 
to be real) while the effective local KS potential v sa {r) 
is the sum of the external, electrostatic and xc terms. 
This potential is modified, for convenience in theoreti- 
cal considerations, by adding an infinitesimal symmetry- 
breaking term to remove any possible degeneracy (cf. 
Ref. H), and by enclosing the system in a large box to 
have fully discrete energy spectrum (and to guarantee 
convergence of spatial integrals). The latter modifica- 
tion of v s(T is not required for bound orbitals and it is 
not used by us in practical calculations of their wave- 
functions and energies with the presently applied numer- 
ical method (Sec. IVII Aft . The spin- up, n-)-(r), and spin- 
down, nj.(r), electron densities are determined with the 
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occupied KS orbitals (p aa (a = 1, . . . ,N a , N-f + N± = N) 
through the relation 

n tr (r) = n[{cf> a(T }](r) = J2^L(r). (1.2) 

a=l 

Since the exchange potential v xa depends only on the 
density n a and it is expressible in terms of cr-subsystem 
characteristics, the dependence on the spin index a in 
V xcn n a , N a , 4>aa etc. can be suppressed hereafter and 
all following discussion refers to one of the spin-up or 
spin-down sub-systems separately. 

The exact exchange potential v x (r) = 5E x [n]/Sn(r), 
defined with exact E x , satisfies an integral equation the 
kernel of which is the static KS linear response function 
involving all (occupied and unoccupied) KS orbitals of 
given spin. This equation was originally derived within 
the optimized-effective-potential (OEP) approach 
and it can be also obtained [ll| when v xa (r) is identi- 
fied as the first-order term in the xc potential expansion 
resulting from adiabatic switching the electron-electron 
interaction in the many-electron Hamiltonian The 
OEP equation for v x can be solved numerically within 
a spatial grid representation for atoms @, H OHl] but 
the involved determination of the KS response function 
matrix makes this approach inefficient for larger systems. 
The solution of the OEP equation with a finite basis for 
orbitals and an auxiliary basis for the response function 
[Til [l7l is usually very troublesome for atoms and 
molecules since it requires careful balancing of the two 
bases to avoid unphysical oscillations of v x (r) f]~Ql — Tsij j . 
Similar problems due to basis imbalance are encountered 
[22l - |25j when the exact v x , expressed in an auxiliary basis, 
is found by following directly the OEP approach, i.e., by 
minimizing the HF-like energy expression for total energy 
given in terms of occupied (KS) orbitals obtained with a 
local potential v s ; see Ref. [26|. An alternative method of 
solving the OEP equation in orbital basis representation, 
which avoids the use of an auxiliary basis by expressing 
potentials in terms of the products of occupied and un- 
occupied KS orbitals, has been proposed very recently 

na. 

A different approach to determination of the exact 
v x , introduced by Kiimmel and Perdew (KP) [28], uses 
the OEP equation rewritten in terms of the occupied 
KS orbitals and the corresponding orbital shifts (OS) 
@, HU; see Sec. UTTXl The latter are solutions of N 
non-homogeneous linear differential equations, hearafter 
called the KP equations, which depend on the unknown 
v x . The so posed problem of determining v x (for a given 
density) has so far been solved in two ways (Sec. IHI CI) : 
(i) by iterating v x and the OS [28[, (ii) with our previ- 
ously developed algorithm [3(| where v xa is found with- 
out any iteration, using the solution of the system of 
equations (algebraic and differential) obtained by a suit- 
able combination of the OEP and KP equations. A fur- 
ther discussion on the exact exchange potential and its 
orbital-dependent approximations, methods of their de- 



termination and applications can be found in the recent 
review [31[ by Kiimmel and Kronik. 

In the present work, we determine the exact exchange 
potential with a novel iterative method based on the rep- 
resentation of approximate v x (at each iteration step) 
in terms of the corresponding OS and energy shifts 
(ES). Such representation was previously derived in Refs. 
I9H28I and [29| for the exact v x , but it is now recognized to 
be an identity relation holding for any local multiplica- 
tive potential (Sec. [XXJ) , in a similar way as the relations 
previously found in Ref. .32!. Our new algorithm, summa- 
rized in Sec. EH consists of (i) calculating the OS and ES 
with a finite-difference formula from solutions of the KS(- 
like) and Hartree-Fock(HF)-like equations (Sec. IIV|) . and 
(ii) iterating v x (for a given density) by modification of 
its ES and OS terms with an appropriate use of the OEP 
equation (Sees. IIII D[ fV)) . (iii) iterating the KS potential 
v s to obtain the convergent density. Thus, solution of 
the KP equations is avoided, while the iteration of v x is 
done in a different way than in Ref. [28]. The method 
for calculating of the OS and ES follows a recent finding 
that the OS are well represented by the differences of the 
respective occupied KS and HF orbitals, see Table I in 
Ref. HH (note that in the latter work HF orbital-specific 
exchange potentials are shown to map very accurately 
onto v x and thus the outstanding proximity between the 
HF and KS orbitals is explained). 

The efficiency of the proposed algorithm is tested (Sec. 
IVII[) for several closed-(sub)shell atoms from Be to Kr 
within the exchange-only KS scheme where the correla- 
tion energy E c is neglected together with the potential 
v c , a part of the xc potential, thus leaving only v x . In 
the performed tests, the obtained exchange potentials, 
which are, in fact, very accurate approximations to the 
exact v x , are compared with the exact exchange poten- 
tial (used as a reference) determined with extremely high 
accuracy using our non- iterative algorithm [30(. Similar 
comparison is done between the corresponding electron 
densities and orbital energies. 



II. IDENTITY SATISFIED BY ANY LOCAL 
POTENTIAL 

We find it useful to begin our presentation of a new al- 
gorithm for the exact exchange potential v x (r) = v^(r) 
by our reformulation of the investigation of Krieger et al. 
Q, continued next by Grabo et al. (29[ and by Kiimmel 
and Perdew [281 ] . It is given in a form of an identity 
satisfied by an arbitrary local (multiplicative) potential 
v x (r) [some restriction on it will be indicated below Eq. 
(|2.9c[) ] . This identity is written in terms of the occupied 
KS orbitals {(/> a (r)}% =1 , Eq. (fTTTaj) [directly and via the 
density n(r), Eq. (I1.2[) ]. auxiliary constants — the en- 
ergy shifts (ES) {-Daal^Li and auxiliary functions — the 
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orbital shifts (OS) {8<fc(r)}% =1 , namely 
where 



(2.1) 



1 N 



JV 

y 

n(r) ^ 

v ; a=l 

is known as the Slater potential, while 



JV 



where 



v* S [{D L aa }A<f>aW) = v* SL (r) = $>r(r)I^ , 

a=l 

(2.3a) 

<t>l{r) 



vTir) = 



n(r) 



(2.3b) 



and 

v° S [{H h a }A<t>a}]{r) = « x OSL (r) = f> x ° a S (r)^( r ) = 

a=l 

1 * 

-7-t E{^( r ) - ^»(r) *»^(r)} , (2.4) 

' ' a=l 

are the ES and OS components of the potential. The OS 
operator v xa s [occurring in Eq. (|2.4[) ] includes a differen- 
tial operator i(r), while the Fock exchange operator 
[occurring in Eq. (|2.2[) ] is a non-local integral one defined 
by 



%[{ct>a}](r)v(r) = -J dV£ 



0a(r)0a(r')v(r') 



r — r' 



(2.5) 

(2.6) 
(2.7) 

where for any potentials difference <5w(r) the ES D aa [5v] 
are defined as its matrix elements 



The ES and OS for the potential are given by 
D% a = D aa [i%-<j*], 
^(r) = ty„[i£-£](r), 



(2- 



while the OS 6<p a [Sv] (r) are solutions of the differential 
equation involving Sv(r) 

£(r)+« s (r)-e a ]^ [a](r) - -w a [Sv](r) 4> a (r) , (2.9a) 

where the right-hand side operator is 

w a [Sv]{r) = Sv(r) - {(f> a \8v\^> a ) . (2.9b) 

These solutions must satisfy the same boundary condi- 
tions as orbitals. In addition, they must satisfy the or- 
thogonality requirement 



The identity (|2.1[) concerns only such potentials t£ that 
the solution 54> a [v^ — v£] of Eq. (12.91) exists. Note that 
(4'a\'Wa\4>a) — 0; an d that u> a [<5u] is invariant with respect 
to a constant shift c of its argument, therefore 



6<j) a [8v + c](r) = 6(/> a [6v](r) . 



(2.9d) 



As we see, all terms on the right-hand side of the iden- 
tity (|2.ip can be evaluated knowing only the occupied 
KS solutions {0 a , e a },?Li and the KS potential v s . The 
dependence on enters via the functionals ES and OS, 
Eqs. (|2.6[) and (|2.7p . The dependence on r enters the sep- 
arate components via {(frail")}, n[r), v^(r) and {S(j>^(r)}. 
The identity (|2.1I) is similar to other three identities ob- 
tained by us in our paper [32J . Any local potential poten- 
tial expressed according to the the identity (|2.ip will be 
named to have the canonical form. This identity will be 
used by us for representing an approximate exchange 
potential in the process of iterative improvement leading 
to the exact one. 

The validity of Eq. (|2.1[) can be proven quite simply. 
For 5v = w x — t) x one should multiply both sides of Eq. 
(|2.9ap by <fi a (i'), sum up over a from 1 to N, use defini- 
tions (|1.2p and (jl.lbp . and utilize Eq. (|l.lal) to replace 
£q <t>a by the equivalent expression. Locality of v x and 
v s is to be taken into account. It should be noted that 
the requirement (|2.9cp is not involved in this procedure. 
However, it defines uniquely the particular solution of Eq. 
(|2.9ap that allows for expressing 6<p% in the equivalent 
perturbation-theory form, Sec. |lVl A general solution is 
5<p^ + c a (j> a , with an arbitrary constant c a . 

There are interesting properties of the canonical rep- 
resentation of a local potential. The replacement <5<Zj„ — > 
8(j)\ + c a (f> a (with c a — a constant) does not change the 
OS potential, Eq. While the OS, Eq. p7T) . are in- 

sensitive to the replacement — > + c x (with c x — a 
constant), see Eq. I|2.9d|) . the ES, Eq. (|2.6p . are changed 
by this replacement 



(2.10) 



Therefore, due to the obvious identity [see Eq. (I2.3bp ] 



JV 



Vr, ^<» = 1, 



the ES potential, Eq. (|2.3ap . transforms as 
, x ES (r)^,f(r)+c x L . 



(2.11) 



(2.12) 



Thus, the identity (|2.ip remains valid for the potential 
u x replaced by (u x + c x ) on both sides of the identity. 

It proves convenient for further applications to fix the 
additive constant of an approximate exchange potential 
v x ( r ) by the replacement 



D , 



D, 



(2.13) 



resulting in 



{4>a\84>a[5v])=Q. 



(2.9c) 



(2.14) 
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With Eq. (I2.14P satisfied, the potential 
v^ s [{D^ a },{cj) a }](r) is exponentially small in the 
large-r region for a general direction r/r (i.e. for r 
not laying on the nodal surface of </>jv(r), if this orbital 
possesses such surface extending to the asymptotic 
region, see Delia Sala and Gorling, [34[ )■ The Slater 
potential v^ l (r) is known to behave asymptotically 
as —1/r. Thus, the asymptotic behavior of the OS 
potential is needed to find the asymptotic form of i£ (r) 
from its canonical representation. 

For the exact (E) exchange potential, = v®, its 
OS part i;° SE (r), Eq. (|2.4I) . exponentially vanishes at 
r — > oo, so it is asymptotically smaller than w|'(r). This 
is a consequence of the OEP equation [see Eq. p.2aj) be- 
low] which implies the OS S<fif f _ 1 and Stfifj in the large-r 
region must behave (with accuracy to their leading ex- 
ponential factors) as —Cn4>n and Cn4>n-i, respectively, 
in order to satisfy Y,a=i 2( t>a = 0. Then u° SE will be 
dominated by the terms with a = N — 1 and a — N, each 
as small as 4>n-i/4>n- 

However, for an arbitrary local potential u L , which 
does not satisfy the OEP equation, the respective OS 
potential w^ SL (r) may be finite at r — > oo. In this 
case, imposing the condition (|2.14[) — though it results 
in the asymptotically vanishing term v E SL (for a general 
direction) — leads to v£(r) that tends to a finite value 



i£(oo) 



,,OSL 



(oo). 



III. 



EXACT EXCHANGE POTENTIAL 



A. Density shift and the OEP equation 

Considering the exchange potential to be a functional 
of the given density of the KS system, we will suppress 
often the dependence of various objects on the fixed KS 
characteristics {(f>b, £b}, v s . As a tool to define the exact 
exchange potential, we introduce the (first-order in OS) 
density shift (DS) induced by <j) a -> 4> a + 54% in Eq. ([j~2")l 

N 

5n L (r) ee 5n[{6<ft}](r) = V 2 a (r) <J#(r) (3.1) 



N 



ee 5h[v^](r) = ^20 Q (r) ty„[t£ - v\\{r) , 

a=l 

where the OS 5<t> a are defined in Eqs. (j2~7|l and (gTSJ). 

As well known, the exact (E) exchange potential, to be 
denoted here v x (r) = (r), can be found as the solution 
of the so called "optimized effective potential" (OEP) 
equation @, H, H 



Vr, 6n[v®](r) = Q. 



(3.2a) 



But this solution can be determined only with accuracy 
to an additive constant. Therefore, to make it unique, 
the additional requirement on v E is imposed 

lim v E (r) = , for a general direction r/r . (3.2b) 



Then D^ N — is satisfied. 



B. Singular components of OS potential 



Since the identity (|2.1| will be helpful in finding the 
solution of the OEP Eq. (I3.2a[) , we make use of the DS 
5n(r) to rewrite the OS potential, Eq. (|2.4p . in two forms: 



1 / N 

= «r) E(- V (^) V ^))} 

^ ^ \ a— 1 



(3.3a) 



v^[m},Sn b ](r) 
1 / N 

= ^) ( Ei^Mr) - (V^(r))-V}^(r) 



(v s (r) + ±i)Sn h (r) 



(3.3b) 



equivalent to the original OS potential only when Sn a 
and Sn h are appropriate [i.e., given by Eq. (|3.ip ] 



OSa 



EE V? Sh 



M}, «n [{$#}]] (r) 
[{6<t>%6n[{8cf> h a }]](r). (3.3c) 



For any u L (r) represented by the identity (|2.1j) . finite 
everywhere, three components vf 1 ^), v^ s [{D^ a }](r), 
v® s {{5(f)^}](r) are also finite when calculated with the 
KS characteristics. But this may be not true for sepa- 
rate terms of v^ S represented in Eq. (|3.3[) . Coulombic 
singularity of the KS potential at each nucleus position 
results in kinks of KS orbitals and density (Kato cusp) , 
and also in kinks of 8(j)\ (r) [when determined from Eqs. 
([2"77]) and (|23)) ] and, therefore, in kinks of 8n[{5(fa}}(r). 
It can be easily verified by expansions (in displacements 
from the point of singularity) that v s {r) 8n[{5(fi^}]{r) 
and i t(r) Sn[{S4%}}(r) taken separately are singular, 
while their sum (v s (r) + \ t(r) ) Sn[{S<j)^}}(r) remains fi- 
nite due to cancelation of singularities. Therefore the 
finite potentials u° Sa and v® sh with arguments indi- 
cated in Eq. (13.3c[) behave differently after truncation by 
means of inserting 5n = 0. The truncated OS potential 
v^ Sa [{S(j)^},0}(r) is singular at nucleus position (because 
the neglected term was singular), while v° sb [{c)0 L }, 0](r) 
is finite (because the neglected term was finite). 



Two known approaches to solve OEP equation 
in terms of orbital shifts 



There are known two approaches to solve the OEP 
Eq. (|3.2al) for ti^(r) by evaluating the OS as auxiliary 
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functions: the iterative procedure, and the one-step pro- 
cedure. In the last case, conceived by us and described 
in Ref. [3(1 one solves the following system of (2N + 1) 
equations: (i) one functional Eq. (|3.2a[) [with Eq. p.l[) 
inserted] , (ii) one functional equation [Eq. (|2.ip for L=E] 

v*(r)=vf(r) + v* s [{D*J}(r) + v$ s [{8<l>*}](r), (3.4) 

(hi) (N-l) functional differential Eqs. ([23a)) . (iv) (iV-1) 
scalar Eqs. (|2.9c[) [in (hi) and (iv), for a = 2,3, ...,N, 
and for L=E], (v) one scalar equation D^ N — 0. [Since 
Eq. (|3.2aj) is included, the OS potential in the form 
v2 Sh [{S(j)a},0}, Eq. (|3.3bj) . is the most convenient.] The 
enumerated system of equations can be solved uniquely 
for one function v^(r), for N functions 5(j)^(r), a = 
1,2,..., N, and for N constants D® a , a = 1,2,..., N. 
These constants satisfy Eq. (|2.6|) . 

The iterative approach was devised, in fact, earlier 
than ours, by Kummel and Perdew (KP) [28]. The im- 
proved approximate exchange potential v" cw (r) is calcu- 
lated from the current approximate exchange potential 
w° ld (r) according to the equation 

< w (r) = «*(r) + v™[{D^}](r) + « x OSa [{^° M },0](r) 

(3.5) 

(with = imposed). This equation follows from 

the identity (|2.ip (valid for arbitrary local potential) 
"spoiled" by inserting the target relation 6n a (r) — 
8n(r) =0 — the OEP Eq. (|3.2aj) for the exact exchange 
potential — into the version v® of the OS potential. 
As follows from our analysis given below Eq. (|3.3c[) , the 
OS term in Eq. (JOJ) calculated for v° ld (r) ^ v®(r) (i.e., 
before convergence) is divergent at nucleus position. KP 
do not comment this feature. According to KP, the it- 
eration cycle based on Eq. (|3.5I) converged quite well to 
the exact potential for extended electronic systems, but 
there were numerical difficulties with evaluating the OS 
potential term for finite systems in the large-r region. 
Therefore alternative recurrence relation was proposed 
and applied by KP 

< ow (r) =v: ld (r)+c6n[{5cf >a [v^-vl}}}(r), (3.6) 

(here a constant c > is a system dependent param- 
eter), see Eq. (13.1[) for the density shift definition. At 
convergence, the second term vanishes, which indicates 
that the obtained potential satisfies the OEP Eq. (|3.2ap . 
In each step of the cycle, the OS are calculated by solving 
a system of differential equations, see Eq. (|2.9I) . 

D. New approach proposed 

Proposed by us in the present paper method of deter- 
mining the exact exchange potential is also iterative and 
makes use of the identity (|2.ip . Three essential novel- 
ties are introduced: (i) calculation of the OS that avoids 
solving the specific differential equation — it will be de- 
scribed in Sec. IIV1 (ii) variational determination of the 



optimal ES described in Sec.[V] and (hi) a new recurrence 
relation for improving potential, where the OS poten- 
tial is evaluated with the help of the modified OS vector 
(Scpi 1 -, <5</>2~S • • • i ^iv") — tne component of the origi- 
nal OS vector {8<p\, 8(fy , . . . , S(fy) perpendicular to the 
orbitals vector (<pi, fa, ■ . . , 4>n), namely, Vr, 

8<fc x (r) = *#(r) - Mr) 8n[{S$}](r)/(2n(r)) . (3.7) 

By construction, the modified OS satisfy the OEP equa- 
tion 

N 

Vr, 5n[{S^ ± }](r)=Y / ^a(r)S4> h a X (r)=0. (3.8) 

a=l 

Therefore, at convergence (L=E) , the second term in Eq. 
(I3.7p . i.e., the modification to 5(f%, vanishes. Of course, 
before the convergence is achieved, S^ 1 - does not satisfy 
Eq. (|2.7p . therefore the identity (12.11) is "spoiled" by using 
Scp^ instead of 54>\. The proposed by us recurrence 
relation 

< ow (r) = vf{r) + v^[{Df a d }](r) + v° s {{5^ d± }}(r) 

(3-9) 

converges fast but only when it is used in combination 
with another step in the recurrence formula — the men- 
tioned optimization of ES; otherwise the convergence is 
very slow or the iterated v x diverges. The relation (|3.9p 
modifying the OS term is a necessary ingredient in our 
approach since the potential u x with optimized ES cannot 
be further improved by mere repeating the ES optimiza- 
tion (see Sec. IV Bp . 

Due to the satisfaction of Eq. (|3.8p . the OS potential 
entering the expression (|3.9p for w x ow can be evaluated 
more conveniently as u° Sb [{^° ld - L } , 0], Eq. (|33b]) . the 
equivalent form, according to Eq. (|3.3cp . There are two 
main advantages of using the modified OS: (i) the OEP 
equation is taken into account in the original (full) OS 
potential expression at each step, and (ii) this modifi- 
cation reduces numerical difficulties (faced by KP) with 
evaluating the OS potential term. 

IV. ORBITAL SHIFTS AND ENERGY SHIFTS 
BY THE NEW METHOD 

It can be easily checked (see, e.g., in 29]) that the 
definition of the OS given in Eq. (|2.9p is equivalent to 

oo 

5<j> a [5v}{r) = M r ) c ja[8v], (4.1a) 

3=1, j¥=a 

where 

c jo [*«]=-«, (4.1b) 

D ja [5v] = {<t>j\5v\<t> a ) . (4.1c) 

It is worth noting that the ES defined in Eq. (|2.8p agree 
with the definition (|4.1c|) of Dj a \8v\. 
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But the expressions given in Eqs. (|4.1[) for {OS , ES}= 
{8<p a [Sv](r),D aa [Sv]} can be interpreted as the first- 
order perturbation-theory corrections to the KS solu- 
tions {cj> a (r) , e a } , stemming from the perturbation Sv(r) 
of the KS Hamiltonian h s (r). So they can be evalu- 
ated also as the first derivatives, with respect to the 
perturbation-strength (real) parameter A, of the so- 
lutions {(f> x [Sv](r) , e^[5v}} of the following one-particle 
Schrodinger equation — the perturbed KS Eq. (|l.la|) : 

{h s (r)+X Sv(r)} ^ [6v] (r) = e x [5v] ^ [5v] (r) , e x < e x a+1 , 

(4.2) 

namely 

6M5vKr)=(dti[6v](r)/d\)\ x=Q = J a [fl)](r), (4.3) 
D aa [Sv]= (de x a [Sv]/d\) | A=Q = °e a [Sv] . (4.4) 

The crucial for our algorithm Eq. (|4.2I) for Sv = — 
will be named the HF-like equation, because its Hamil- 
tonian can be viewed as a modified HF Hamiltonian: 

h s (r)+X(v^(r)-vl(r)) 

= t(r) + (v s (r) + Ai£(r)) - \v*(r) , (4.5) 

containing, besides a local potential, the (scaled by 
—A) characteristic Fock exchange potential operator, Eq. 
(|2.5[) . Therefore, after minor adaptation, standard HF 
codes are suitable for numerical solving Eq. (|4.2I) at any 
finite A. We will also calculate the OS and ES generated 
by some local perturbing potential 5v(r) instead of non- 
local Sv(r). The corresponding perturbed KS equation 
[like Eq. (I4.2[) . but with Sv replaced by Sv] will be named 
the KS-like equation. 

Having solutions of Eq. (|4.2p for a series of A values, 
e.g., A = . . . , —2k, — k, 0, k, 2k, . . . , (the A = solutions 
mean the original KS solutions) the derivatives of orbital 
energies, Eq. (I4.4p . can be found with sufficiently high 
accuracy by applying a finite-difference (FD) formula of 
numerical differentiation: 

e a [Sv] = ^e a [5v] + 0(K p ), (4.6a) 

M°e a = (e«-e a )/ K , (4.6b) 

M e a = (e:-e- K )/(2K), (4.6c) 

^e a = (-2e- K - 3e a + 6e* - ef )/(6/s), (4.6d) 

^°e a = (8e«- 8e- K - + e-^)/{l2n), (4.6e) 

and so on. The derivatives of orbitals, Eq. (|4.3[) , can be 

o 

found from analogous expressions for a particular ' p,K V a 

kM(r)= ^i a [6v](r) + 0( K n- (4-7) 

Here p represents the number of the HF-like equations 
to be solved, and, simultaneously, the power of k in the 
error term, see Eqs. (|4.6a|) and (|4.T[) . 

When performing numerical differentiation of (j>a( r ) 
with respect to A, one must take care to choose con- 
sistently phase factors of involved wave functions with 



various A. If </> A is a real eigenfunction of Eq. (14.21) . 
then (—4>a) is also a valid real eigenfunction. To see 
if the given </> A is close to <j>°, let us express it as fia — 
C (4>a+S4>a) (the allowed C is +1 or — 1) and calculate the 
integral 7 A = / d 3 r {(j) X -(jP a ) 2 = 2(1-C)+C J d 3 r(^) 2 
(</>° and </> A were assumed normalized). We see that 
when <fia and 0° match (C = +1), the result is I x = 
J d r(S(j)\) 2 , while in the opposite case (C = —1), we 
find I x = 4 — J d 3 r (<50 A ) 2 . Taking a 'conservative' es- 
timate that < Jd 3 r(<5^) 2 < 2, we recommend to 
transform </> A — > —(j) X if I x is found to exceed 2. Then 
and </>° can be partners in FD formulas. 



V. VARIATIONAL DETERMINATION OF 
MODEL POTENTIAL 
PARAMETERS 

A. Minimization of density-shift-finiteness 
indicator 

While the exact exchange potential v^(r) is defined 
by vanishing of the density shift in the whole space, Eq. 
(13.2a[) . the density shift due to some approximate ex- 
change potential v^(r), Eq. (|3.1[) , remains finite. As 
a convenient global indicator of density shift finiteness 
(DSF), we define the functional 

A D sf[<H = || MU ) (S- 1 ) 
in terms of the functional norm 




The weighting function w(r) > may be helpful in select- 
ing regions of special importance for the DSF indicator, 
e.g., w(r) = n p (r)/ J d 3 r' n p (r') enhances contributions 
in the low-density regions when 1 > p > 0. We use 
w = 1/N 2 (unless stated otherwise). 

Let us construct a model (m) exchange potential 
as a function of some parameters. Since, obviously, 

A DSF [Sh[vf]] > A DSF [Sn[v*]] = , (5.3) 

the minimization of Adsf with respect to these 

parameters leads to an improved potential, closer to , 
because the corresponding Adsf becomes closer to zero. 

B. Optimization of TV 1 parameters 

Having some approximate exchange potential writ- 
ten in its canonical form, Eqs. (|2.1[) - (|2.4[) . we are going 
to define a model exchange potential u™ by replacing the 
^-dependent ES {D^ a }% =1 with the model ones {D™ a } 
(they will play the role of variational parameters for the 
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mentioned minimization; at convergence, they will be- 
come {D*,}): 

N 

<(r) = v s x \r) + ]T vT(r) ^ + « x ° SL (r-) . (5.4) 

a=l 

In fact, can be expressed as a corrected v^, namely 
<W=^(r) + ^< a (r)Ca (5-5) 

a=l 

where {d,^} can play the role of new model parameters, 
replacing previous ones by 

C = ^o " Daa ■ (5-6) 

The task of finding the model ES {D™ a } that minimize 
the DSF indicator A D SF[^ m ] for 6n m = 5n[{S4> a [v^ - 
fx]}] j Eq. p. II) . is equivalent to minimization of this 
Adsf with respect to parameters {d™}, Eq. (|5.6p . since 
the ES {-Dq Q } are fixed during the minimization. 

Due to the linearity of the left-hand side of the dif- 
ferential Eq. (|2.9ap in the OS and the linearity of its 
right-hand side in the potential shift, we can evaluate 
the total OS [yielding 8n m from Eq. (|3.1[) ] as a sum of 
their constituents 

JV 

6<t> a K ~ €}(r) = 6<fc(r) + £ 5^{r) d£ , (5.7a) 

b=l 

where the partial OS, i.e., 

<ty£(r)=ty«[t£-t£](r), (5.7b) 

Scb b a b (r)=SM< b ](r), (5.7c) 

are the solutions of the differential equation (|2.9a|) in 
which the functional argument of w a [•] is the same as the 
argument of Scj) a [.] seen in Eq. (|5.7b|) or ()5.7c[) . In prac- 
tice, the partial OS are determined from the solutions of 
the perturbed KS equations, see Sec. IIVI 

The optimal parameters {d™ = d^ a }, determined by 
minimization, are solutions of the set of algebraic equa- 
tions dA DSF [5n}/dd™ = 0, i.e., 

N 

^ ac d c L c = «£, a=l,2,...,JV, (5.8a) 

c=l 

where 

A ac = J d 3 rw(r)6n aa (r)6n cc (r)=A ca , (5.8b) 

*~/d-r«(,)*FW*Kr), (5.8c) 
with [cf. Eq. ([3T]) ] 

N 

( 5n L (r)=^20 a (r)^ (r)j (5 . 8d ) 

o=l 
N 

6n cc (r)=J2^a(r)SVa C (r). (5.8e) 

0=1 



Let us note that E&=i = 0. The following steps lead 
to this result: £ 6 *C = £ 6 ty.[«?] = ^.Erf] = 
<50 a [l] = 0, due to identity (pOTJ) : therefore £ 6 <5n bb = 
jC a ^^a b = 0- With the singular coefficient ma- 

trix, det({A oc }) = 0, the solution of the system (|5.8a[) is 
not unique. By choosing (arbitrarily) dj} N = and tak- 
ing a = 1, 2, ... , (AT — 1), we obtain a truncated system 

Ad L = b L (5.9) 

having a unique solution 

d L =A _1 b L . (5.10) 

Here A is (JV - 1) x (JV - 1) square matrix {Ac}^Tii, 
while b L and d L are (JV— 1) x 1 column matrices {b^}^^ 1 
and {d^ c Si. The possibility of arbitrary choice for 
d NN reflects the fact that optimized v™(r) can be found 
only with accuracy to an additive constant. 

Finally, the approximate exchange potential opti- 
mized with respect to (JV — 1) parameters is 

Af-l 

v^ N ~ l (r) = ^(r) + £ d a L a . (5.11) 

a=l 

Because of the factual definition (|5.5p of the model poten- 
tial, no particular (like canonical) form of was needed 
to obtain the last result. At the given step of iterations, 
we need to calculate only the column b L (as dependent 
on i£), while the matrices A and A -1 can be prepared 
in advance (as common to all iterations at fixed den- 
sity). From the above form of w^'^ -1 it follows that 
[u x J,Ar_1 (r) — i>x(r)] i s exponentially small in the large-r 
region for a general direction r/r. 

Once the optimized model exchange potential v^' 1 ^^ 1 
is found, one may be tempted to repeat immediately such 
optimization procedure by starting with i;^'^ -1 as a po- 
tential playing the role of v^. However, the correspond- 
ing new model potential u™ with the new variational 
parameters d™ could then be represented, according to 
Eq. ([53]) with dj} N = 0, as 

N 

<\r)=vt N - 1 (r) + J2< a (r)d^ 

a=l 

N 

= u£(r) + £ <»« L Q + Co') ■ (5-12) 

a=l 

Since the above potential has exactly the same func- 
tional form as v™, Eq. (|5.5p . but with parameters d™ 
replaced by (d^a + d™ ), the indicators Adsf correspond- 
ing to and to attain the same minimum, first one 
at d Z = d Li second one at d£ = d£ = 0. This means 
that the two optimized model potentials are identical, 
iij ■ N ^ 1 = v^-^^ 1 , so no further improvement of u™ can 
be obtained by repeating the described minimization of 
Adsf- 
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C. Optimization of N parameters 

As an improvement to the model potential (|5.4[) used in 
Sec. IV Bl we propose to include a new parameter for op- 
timization — the "strength" C OSm of the model OS po- 
tential, in addition to N parameters — the ES {-D^j^Li- 
Namely we define a model exchange potential 

N 

v?(r) = v s x \r) + £ < a (r) £>™ + C OSl X° SL (r) (5.13) 

a=l 

by replacing in Eq. (|2.1[) the L-dependent objects 
{Pt)a=v v? Sh (r), with the model ones {D^}^, 
C v® (r), which, at convergence, will become 
{D^}, v x SE (r). It will be convenient to change the no- 
tation 

D™^C os ™, v™(r)=v° Sh (r), (5.14) 
to have the model exchange potential in a simple form 

N 

S?(r) = vf(r) + vT(r) K a . (5.15) 

a=0 

which can be further rewritten as 

N 

v™(r) = w£(r) + J2 <» C • (5.16) 

where d™ a are given analogously as in Eq. (|5.6[) , now valid 
also for a — if we set D^ Q = 1. The best model param- 
eters are determined by minimization of the DSF indi- 
cator Apgp[5n], The same equations as in Sec. IV Bl are 
obtained, only the range of indices is extended now from 
the initial a,c — 1 to a, c = 0. Due to singularity of 
the coefficient matrix, we again choose dj} N = 0. The 
following (truncated) system of equations 

A' d" L = b' L (5.17) 

is obtained, having a unique solution 

d' L = (a')~V l . (5.18) 

Here A is N x N square matrix {A ac }^ ~2 , while b L and 

d' L are N x 1 column matrices {b^}^ 1 and {c^}^ 1 . 
It should be noted that the elements Aq c — A c q and Aqq 
are specific to v^, because they depend on v®° = v® Sh 
via Sn 00 . Therefore, they, b L , and (A ) need to be 
calculated anew at each step of iterations. However, as 
shown in Appendix, costly evaluation of (A ) at each 
step can be avoided when a special method of solving Eq. 
(|5.17j) is applied, using A -1 prepared in advance. 

Finally, the approximate exchange potential v x opti- 
mized with respect to N parameters is 

N-l 

t£» = « x L (r) + J2 <» tia • (5-19) 

a=0 



The OS potential u° SL of the canonical form of v x is 
needed to obtain and write down this result, since it is 
involved in Eq. (|5.14|) . It should be noted that costly 
evaluation of v x L according to the definition (|2.4I) can 
be avoided, if the identity character of Eq. (|2.1I) is used. 
Then 

N 

v° Sh (r) = i£(r) - v s x \r) - ]T <» D\ a . (5.20) 

a=l 

The ES {D^ a } are easily evaluated as matrix elements, 
Eq. flU} with (|2~8)) , or with the help of the method de- 
scribed in Sec. ITV|) . 



VI. ITERATIVE ALGORITHM LEADING TO 
THE EXACT EXCHANGE 
POTENTIAL AND SELFCONSISTENT DENSITY 

A. Iteration improving density: step I = 

When demonstrating efficiency of our algorithm for de- 
termination of the exact exchange potential v E of DFT, 
two loops of iterations leading to self-consistency are 
employed: the external loop for improving the density 
(labeled by {£}, I — 0,1,...) and the internal loop 
for improving the exchange potential (labeled by [k], 
k = 0, 1, . . .). At fixed {£}, all iterations are performed 
for the given density n™' of the KS system (and, there- 
fore, for its fixed characteristics {4>t, ej}, v s ). 

As a starting approximation to the exact exchange po- 
tential at I = we choose the Krieger-Li-Iafrate (KLI) 
potential Q which can be defined as having the canonical 
form with the OS component equal zero 

N-l 

^ LI [n](r) = v s x \r) + £ ««(r) (6.1) 

o=l 

and the ES of the ES component satisfying 

D™ = D aa [v™ - vl] , (6.2) 

[see Eqs. flU} and (JHHJ)] ■ Therefore {D™} can be 
determined from a system of algebraic equations, a = 
V 1. 

N-l 

= D aa K S1 - vl] + J2 K b ] . (6.3) 

b=l 

Note that = (<f> N \v^ hl - vl\<f) N ) = is satisfied. 

The KS solutions are determined next self-consistently 
with this v^ Ll [n] playing the role of the exchange contri- 
bution to the total KS potential v s [n] — w cxt + v cs [2n] + 
v x [n] (the exchange-only approximation of DFT is ap- 
plied by us to spin-compensated systems where n = n-j- = 
ni). In this way for I — the starting exchange potential 
and KS system are constructed with 

v n [n {0}] = v KLi [n{ o }] ( w {0} = ^ + Wcs [2n{°}] + 4°1 , 

(6.4) 
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based on the self-consistent solutions. 

Before entering calculations at (£ = 0, k — 0) one 
should choose values of some general parameters: p, n — 
for numerical differentiation in Eqs. (|4.6[) , (I4.7[) . /co > 1 ~~ 
the final iteration step, a m i x — for mixing KS potentials 
in Eq. (|6.25[) (applied by us value is 0.4), r^DSF and ?7dv 
— thresholds of accuracy [see Eq. (|6.26j) ]. w(r) — weight 
function used in Eq. (|5.2|) and elsewhere; w = 1/N 2 is 
applied here. 

B. Iteration improving exchange potential: 
step k = 

Quantities depending only on the fixed characteristics 
of the actual £th KS system are calculated: (r), Eq. 

w»(r), 50(r) = 5<p a [vf]{r) = M0>f](r), 
Eqs. (I23bl . d4~71) . for a= I,..., N,b = I,..., N -I (KS- 
like equations are to be solved, see Sec. IIV|) : (5n bb (r), 
A 6c , Eqs. dEHeJl, (|5^8b| . for b, c = 1, . . . , N - 1; A" 1 
as reciprocal to A — the matrix {Ab c } b In fact, all 
these calculations can be performed at I = only, [see 
the comment below Eq. (|6.13[) ]. while at £ > solely 
v^(r) and {v^(r)} are to be calculated. Set k = 0. 



C. Iteration improving exchange potential: 
step k + 1 

1 . Initialization 

For given v£\ using Eqs. ((43J), ([iTf]) . (TO)) and (|4~6|) . 
OS and ES are calculated (for a = 1, . . . , N) 

60 = 84> a [vM - <%] = M ki0 ~ C] , (6-5) 

dm = Daa [vW - <] = M° €a[v m _ C] , (6 . 6) 

- the HF-like equations are to be solved, see Sec. IIVI 
Such situation occurs at k = 0. However, for k > 0, when 
the potential is known to be the sum 

v W=v^+6vW (6.7) 

and the quantities 60 1 ' = 5(f) a [vx ^ — v^} and 
D^aa 1 ' = D aa [v^ ^ — v^] are available from the previous 
step, then 

60 = 80~ 1] + M 0^ ] ] > (6-8) 

= Dt~ 1] + M °e a [6v^} , (6.9) 

so the less expensive KS-like equations are to be solved 
because the potential 5v\^ is local. 

The obtained OS and ES are used in further calcu- 
lations of the step. In particular, the ES potential is 
determined as 

N 

vf^(r)=J2<ir)D^, (6.10) 

a=l 



allowing to determine the OS potential directly, see Eq. 
(|OD|) . as 

v° s W(r) = 4 fc l(r) - vf(r) - v^(r) . (6.11) 

In this way all terms of the canonical form of (r) arc 
available. 

Next, a modification of (r) is to be calculated. If k 
is even, the modified ES term is determined: either via 
(N — l)-parameter optimization performed in Sec. lVI C 21 
or via iV-parameter one in Sec. lVI QUI Otherwise (for odd 
k) the OS term is modified in Sec. IVI C 41 

2. Optimization of ES term: (N — 1) -parameter version 

Optimization step discussed in Sec. IV Bl is imple- 
mented. The OS {6 ] } calculated in Sec. IVI C II are 
used. Applying Eqs. fl5j3HJ) , (l5lfc]) . (I5TU1) for L = [k] cal- 
culate 6nW, {b l a ] } : and finally (note that {60} and 
A -1 were obtained at the step k = 0). The optimized 
potential, see Eq. (|5.1ip . 

v™' N - l (r) = 4 fel (r) + £ <°(r)4a ] , (6-12) 

o=l 

represents the [k + l]th potential, which is characterized 
by the modifying term 

N-l 

54 fc+11 (r) = £< a (r)dW (6.13) 

a=l 

(here k + 1 is odd), see Eq. (jB~2Tj) . 

It is worth noting in passing that it is sufficient to 
use 60, 6n bb and A -1 obtained at the starting step 
I = 0, k = (instead of calculating these objects at the 
actual £, k = 0), because this simplification introduces to 
dl fc ] an error Sd^^ linear in a small KS potential modi- 
fication 6vs e \r) = Vs(r) — vj°^(r), resulting, therefore, 
in an error of the DSF indicator minimum Adsf that is 
quadratic in 6d^ k ^ £ \ so quadratic in <5«i^. We observe 
no appreciable effect of this simplification on the conver- 
gence rate. 

Go to Sec. IVI C~5l — the closing subsection. 

3. Optimization of ES term: N -parameter version 

Alternatively to previous Sec. IVI C 21 the optimization 
step discussed in Sec. IV Cl is implemented here. The OS 
{<5</>L fc '} and OS potential v$ S ' fc ', calculated in Sec. lVICTl 
are used. For L = [k] calculate v°°, Eq. (|5.14|) : {60 = 

60 a [v™} = M i a K°}}^ =1 , Eqs. 631), 6U); 6n°°, Eq. 
(I5.8cl) : Aq c = A c q, Aqq, Eq. (HHbJ); fcff 1 , Eq. (1TO . 
Solve the system of equations 

A'd'M = b'W (6.14) 
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for unknown vector d [see definitions of A , b ^ , d 
below Eq. (|5 . 18[) ] . Inexpensive method of solving this 
system is shown in Appendix, where the matrix A -1 (cal- 
culated once at the k = step) is used. 
The optimized potential, see Eq. (|5.19[) . 

4 fe] '» = 4 fel (r) + £ «£» d «« > ( 6 - 15 ) 

represents the [k + l]th potential, which is characterized 
by the modifying term 

N-l 

Sv^ +1 \r)='£<(r)di k J (6.16) 

a=0 

(here k + 1 is odd). The simplification mentioned below 
Eq. (I6.13|) concerning 8<jr a , 8n bb and A -1 is applicable 
here too. 

Warning: at the starting step (£ = 0, k = 0) the 
iV-parameter optimization is impossible because i>°° = 
ijOSKli = 0. This is a specific property of the chosen 

initial exchange potential v$ = v^ Ll , see Eq. (|6.1|) . The 
(N — l)-parameter optimization should be applied in- 
stead. 

Go to Sec. IVI C~5l — the closing subsection. 

4- Modification of OS term 

The recurrence relation, Eq. p.9[) , proposed earlier in 

Sec. MI D| is implemented. We transform the OS Scj)^ 
[obtained in Sec. IVI C 1| into the modified OS according 
to Eq. ([37f)l , Vr, for a = 1, . . . ,N: 

S^ x (r) = 5c^\r)-<f> a (r)8n^{r)/(2n^(r)) , (6.17) 

to obtain the new OS potential 

„os[*+i] (p) =v™ h [{8ct>W ± },0](r) (6.18) 

[Eq. (j3.3bl) ]. for the use in the recurrence relation 

(r) = vf(r)+vf[{D^}](r)+v^ k ^(r) . (6.19) 

The modifying term is therefore 

Sv^ (r) = v° s ^ (r) - v° s ^ ( r ) ( 6 .20) 

(here k + 1 is even). Note that the OS potential i^? S ' fe ' 
was calculated in Sec. IVI C II 

5. Closing 

The new exchange potential is obtained as 

„[*+!] ( r ) = „[*] (r) + Sv [k+i] (r) _ ( 6 21) 



The steps for optimization of the ES term and modi- 
fication of the OS term have implemented fully our new 
ideas how to update the approximate exchange potential 
vH^ (r). Therefore, these steps should be repeated to gain 
further improvement. Due to the modification of the OS 
potential, optimization of the model ES parameters can 
be effective in the next k step. 

Perform the replacement k :— k + 1 to update the 
step number. If k ^ kg, go to Sec. IVI CI to perform a 
consecutive k step, otherwise continue as below. 



D. Iteration improving density: step £+ 1, 
termination 

When the exchange potential iterations are terminated 
at k = ko (quite small values can be taken, like ko = 2 

or fco = 3) the accuracy of may be still low, but it is 
already sufficient to improve the density in the next, (£ + 
l)th step. This accuracy can be measured by (smallness 
of) 

A| F EA DSF [W fc '] (6.22) 

[see Eq. (|5.ip ], where Sn^ is calculated from {50^} (ob- 
tained in Sec. IVI C II) by applying Eq. (|5.8d[) for L = [k]. 

In order to fix the arbitrary additive constant of the 
obtained w x fc °' in agreement with Eq. (|2.13[) . calculate [see 
Eqs. (HU) and dUHJl] 

^ = <^I4*° 1 -*II^>, (6.23) 

and, next, noting that X)^=i w x°( r ) = 1> perform the 
replacement 

4 fc »](r):=4 fc «l(r)-^. (6.24) 

For the (£ + l)th step, the potential v^ o] [nW](r) is 
used to construct the new KS potential as a mixture of 
potentials (regulated by a parameter < a m i x < 1) 

= a miK (v CKt + v„[2nW] + vf ] [n«>]) + 
(1 - a mix ) . (6.25) 

The KS solutions with this KS potential are obtained. 
They define the density n^ +1 ^(r). 

At this point one should check, using some criteria, if 
the density n^ e+1 ^ approximates the selfconsistent den- 
sity within some presumed accuracy limits. In our ex- 
ample for this check (discussed later on), the indicator 
of density values (DV) discrepancy Adv^i,"^] = — 
n 2 \\ w is introduced and A^+ 1}W = A D v[n {e+1} , »W] is 
calculated. If 

(Ag» ] F < ?7dsf) A (A^+ 1}W < mv ) (6.26) 

is true, iterations over £ are to be terminated. Here the 
parameters < ryDSF <C 1 and < 7/dv <C 1 represent 
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the chosen thresholds of accuracy. The determined ex- 
act exchange potential is v x (r) w vt' ] {r) for all r [the 
correct large-r behavior is guaranteed due to the replace- 
ment 

But if the expression (I0.26[) is false, iterations over I 
are to be continued. Define the initial exchange potential 
as 

4°1 [n^ +1 >] (r ) = 4 fe °l [#] (r ) , (6.27) 

perform the replacement £ :— i + 1, and go to Sec. IVI Bl 
— the exchange potential iteration step k = 0. 



VII. RESULTS OF TESTS 
A. Numerical method 

The algorithm described in the previous section is 
successfully applied to several spin-compensated closed- 
(sub)shell atoms, from Be up to Kr. Similarly as in our 
previous work [3fjj, the calculations are performed using 
the highly accurate pseudospectral (PS) method in which 
KS and HF orbitals 4> a {r) — R n i(r)Yi m (f) (where a = 
(nlm) is specified with the main n, angular momentum I 
and magnetic m quantum numbers) are represented with 
the values of the radial functions Xnl(r) — r R n i{r) at 
M+2 discrete points (nodes) r*j = r(Xi) which correspond 
to an increasing sequence Xo = — 1, x-y, . . . , Xm, Xm+x = 1 
via the scaling r(x) = L(l + x)/(l — x) (the parameter 
L = 0.3 is used). Simultaneously, the KS- and HF-like 
equations are transformed into eigenvalue algebraic equa- 
tions for the orbital energies e n i and the corresponding 
eigenvectors [uf l ] = [g(ri)xni( r i)]fLi which are solved us- 
ing standard techniques; note that Xni{fi)=0 for tq = 
and tm+i = oo. The form of orbital-independent func- 
tion g(r) and other details on the application of the PS 
method to the KS equation are given in Refs. HH andl36l. 
With the presently chosen scaling r(x), it is also possible 
to fully implement the PS representation in numerical 
solution of the HF(-like) equation for atoms, i.e., to use 
— like for the KS equation — solely the radial function 
values {Xnz(fi)}f=i at tne nodes. This novel technique, 
to be reported elsewhere (39|, differs from Friesner's ap- 
proach [33, [38[ where the HF equation is solved by com- 
bined use of orbital basis and the PS representation of 
applied functions. In the present calculations applying 
the PS method, a very moderate number of M — 120 
nodes is used and it is found to be sufficient for the accu- 
racy 10 -10 Hartree of the occupied orbital energies ob- 
tained for given one-electron Hamiltonian (it is true for 
both the KS(-like) and HF-like equations) while the self- 
consistent orbital energy values are two or more orders 
less accurate due to inaccuracies of v x of the FD origin 
(even at convergence) as it is discussed later on. 

In the PS method the exchange potential v x (r) needs 
to be known only at the nodes r = fj (i = 1,... , M). 
Thus, it is also true for the constituent terms of i> x (r) 



in the canonical representation and to its increments 
determined in the iteration algorithm. The ES term 
v x S ( r i) i s readily calculated with {Xni( r *)}> nl <= occ > 
while the Slater potential wf 1 ( r i) is obtained with all val- 
ues {Xnl(fj)} when expressing the action of the Fock 
exchange operator on the orbitals in the PS represen- 
tation [39j | . Also the radial OS functions 5xni{r) and 
their derivatives d[5xnl(i r )]/dT, used to modify the OS 
term of the iterated u x (r) (Sec. IVI C 41) . need to be de- 
termined only at r = r,. This is done with the (modi- 
fied) FD formula (|4.7|) using [instead of (f>^ (r,-)] the values 
Xni( r i) which are found directly by diagonalizing the ma- 
trix of the HF-like hamiltonian (|4.5[) . In calculation of the 
Slater, ES and OS terms of v x (r), the KS functions xm(r) 
— which are not accurately represented for large r in the 
PS method [3C| — are approximated for r > roi with 
the asymptotic form XniH r ) = C n ir l l linl e~^ nir where 
fini = y— ZtnV, the parameters C n i are found by match- 
ing Xni{r) and xiT^( r ) a ^ ^ ne ra dius r = roi, chosen by 
our code individually for each orbital (nl) (e.g., for Ne, 
roi = 2.85, 11.76, 16.89 a.u. for (nl) = Is, 2s, 2p). The 
effect of similar inaccuracies of the OS determined with 
the FD formula at large r can be avoided by setting the 
exponentially decaying OS term v® s (r) to for r larger 
than some cut-off radius 7-02 (e.g., 9, 12, 14 a.u. for Ne, 
Ar and Zn), without affecting the hi gh accuracy of the 
calculated orbitals and their energies [301 ] . 

The PS nodes rj are also used in (a variant of) the 
highly accurate Gauss-Legendre quadrature applied to 
calculate integrals representing various quantities applied 
in our method, including the DSF and DV indicators or 
the elements of the matrix A and the vector b L used in 
the ES optimization. 



B. Efficiency of present algorithm for exact 
exchange potential 

To test one of the key elements of the present al- 
gorithm, the approximate ES 5t n i = D n i, n i — ' p ' K Wi 

o 

and radial OS 5R nl = M R n i are calculated for the 
exchange-only OEP density n OEP and the correspond- 
ing exact exchange potential i> E [n OEP ] (as well as the 
respective orbitals and their energies) by using the FD 
formulas (|4.6[) and (|4.7[) with different p and k. The ob- 
tained results are compared — in Table U and Fig. [T] 
- with the exact ES and OS, i.e., <5e E z = D°, E ^ and 
<5i?^ ; EP determined with our method of Ref. [3(| For given 
order p of the FD formula (|4.6[) . the approximate ES 

Se nl =M ° £nl [ v v - i,?} and OS SR nl = MR nl [v* - v^} 
become closer to (5e E z and <5i? E z when k becomes smaller. 
For p = 1 and p = 2, the discrepancies Se n i — Sef l{ 
well obey the predicted power law k p , but they de- 
crease so quickly with increasing p for each considered 
k = 0.05, 0.1, 0.2 (Table U) that for p = 3, p = 4 they 
approach the accuracy limit, and, in consequence, their 
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k p dependence is not clearly seen for these p. The accu- 
racy of the calculated ES Se n i (and their accuracy limit) 
is a combined effect of the n p non-numerical error of the 
corresponding FD formula (|4.6p and numerical inaccu- 



racies of leading to an additional error in the ES, 
proportional to 1/k (which is a factor present in FD for- 
mula (I4.6p ). Similar conclusions hold for the accuracy of 
the OS calculated with the FD formula (I4.7[) . Hereafter, 
the OS and ES (as well as their increments) used in de- 
termination of Vk are normally calculated with the FD 
formulas (I4.6[) . (14.71) for p = 2, n = 0.05 (invoking only 
two terms or R^i, with A = — k, k), unless the depen- 
dence of the result accuracy on p and k is studied, as in 
Tables U HID and Fig. Q] Note that the FD formulas of 
higher order p (3 or 4) — though requiring calculation of 
more terms and — should be applied when inac- 
curacies of these terms become larger since this adverse 
effect can be compensated with use of larger k without 
losing the accuracy of the resultant ES and OS. 

The efficiency of the algorithm for self-consistent deter- 
mination of n and v x has been tested for closed- (sub) shell 
atoms from Be to Kr and the selected results for the Ne, 
Ar and Zn atoms are presented in Figs. [5HS] and Tables 
HI! Mil The quality of the approximate exchange poten- 
tial wlf' obtained for a given electron density (n OEP is 
used in this test) is first assessed for the Ne atom by di- 
\k\ 

rect comparison of w x with the exact exchange potential 
«x [n OEP ] (its plot is given in Fig. 1 of Ref. 28(a) ) ■ The it- 
erated potential v$ converges quickly to v E (Fig. Hfc)), 
so that (the maximum magnitude of) the discrepancy 



which is a suitable measure of how well the OEP equa- 
tion (|3.2ap is satisfied by approximating f E . In tests 

for Ne, Ar, and Zn (Fig. [3]), the indicator Aq SF calcu- 
lated for a fixed density (n OEP is used) decreases rapidly 
during iteration, roughly as cfc -7 (c, 7 — constants), 
until a minimum level Ap SF is reached after some num- 
ber k± of iteration steps. This corresponds to obtaining 

convergence of for k = fci since Ap SF changes in 
oscillatory manner for k > k±, with amplitude of the 

rfci 

order of Ag SF . The decrease of A^gp is faster for the N- 
parameter ES optimization than for its (N— l)-parameter 
variant though the very similar Aggp level is achieved in 
both cases. Its value depends on the accuracy of the ES 
and OS obtained with the FD formulas and it is around 
10~ 7 for presently used p — 2, n — 0.05 but it can be 
reduced below 10~ 8 for suitably chosen p, k (e.g., p = 2, 
k = 0.02 or p — 4, k — 0.1). The number of steps 



required to obtain A 



[k] 

DSF 



10 7 , corresponding to con- 
vergent v* 1 in the present test, is k\ — 30, 17, 21 with 
(N — l)-parameter optimization and k± = 13, 13, 11 with 
./V-parameter optimization for the Ne, Ar, and Zn atoms, 
respectively. However, note that the found difference in 
the speed of convergence of — especially evident, in 

our test, for Ne and Zn — occurs only when we consider 

rfci 

Vx calculated for a given density. This advantage of N- 
parameter ES optimization is not retained in the practi- 
cal calculations when both the electron density and the 
exchange potential are found iteratively; see later on. 



A4 fcl 



,[fc] 



is reduced 10, 10 2 , 10 3 and 10 4 times 



after k = 9, 16, 22, 29 iteration steps, respectively, in 



[o] 



i> E (equal to —0.25 Hartree 



comparison to Av : 

for r = 0.26 a.u.). Largest discrepancies | Av^ 1 1 are found 
in the interval 0.2 a.u. < r < 0.35 a.u. corresponding to 
the KL intcrshell region where the characteristic bump 
in the r-dependence of w E is located in the Ne atom. In 

this region the iterated potentials v{ k] change most sig- are f ten i arger t han A^ 1 , despite the mentioned over- 



Note also that for all even k the values Aq SF are 
smaller than A^gp 1 ' (see Fig. [3]) since the correspond- 
ing v\^ are obtained by minimizing Adsf of v* ^ by 
treating the ES (and D 00 ) of this potential as N — 1 (or 



N) variational parameters. For odd k, we find that A' fc ' 



DSF 



nificantly at even iteration steps k when is obtained 
by modifying the OS term of v x k — this leads the sign 
change of A«x in comparison to Ai>x 



The (N - 1)- 

parameter ES optimization of v x k ^ — applied in this 
test for Ne to calculate v x at odd steps k — has much 
smaller effect in the intershell region. However, it leads 
to much larger changes of the exchange potential v^ for 
< r < 0.2 a.u. (the K shell region) which may be com- 
parable to those induced by the OS modification at even 
iteration steps. The figure [5Ja),(b) demonstrates that di- 
minishing of extremal |<W fc l | is in accord with diminishing 

of extremal |Au x fc '|, Fig. [2Jc). 

In practical calculations, when the exact potential v® is 
not known beforehand and thus cannot be used as a refer- 
ence, the assessment of the v^ quality (suggested by the 
comparison of Fig. [5Jc) with Fig. [H[a),(b)) is done with 



aid of the DSF indicator A 



[k] 

DSF 



l<W fe i||, Eq. derm 



all rapid decrease of Aq SF 



If the ES optimization is not used, the corresponding 
DSF indicator (marked with triangles in Fig. [3]) does not 
decrease with k (or decrease very slowly, e.g., for Zn) so 

that the iteration of v x does not converge. In particular, 

rfci 

it is found that, with increasing k, the potential v x oscil- 
lates around u E in the KL intershell region, changing the 
sign of Av^ after every the OS modification step, but 
\Av* : I does not decrease during iteration. This adverse 
behavior cannot be remedied with (often used in similar 
cases) mixing of the current v x with the one obtained in 
the previous step: thus the inclusion of the ES optimiza- 
tion step becomes essential to obtain convergence of the 
iterated exchange potential. 
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C. Efficiency of present algorithm for 
self-consistent determination of density 

Further calculations are performed using the itera- 
tions described in Sec. IVI1 within the exchange-only KS 
scheme. The convergence of the iterated density n — 
is monitored using the DV indicator A D y'^ ^ = 
||tiOT — The discrepancy of from the exact 

OEP density n OEP is measured with another DV indica- 
tor 



A W.OEP = || n{i} _ n OEP||_ j t is fmmd m Qur tegt 



Ne (with fco = 2 and the (N — l)-parameter ES optimiza- 
tion) choosing ?7dsf = i)dv = 10~ 4 , 10" 5 , 1(T 6 , 10~ 7 
leads to termination after, respectively, I — 6, 10, 15, 20 
steps of the density iteration while the resultant be- 
comes increasingly accurate since we find A^y OEP w 

The occupied orbital energies e„/, obtained with the 
KS potential Vs at consecutive iteration steps i 7 con- 
verge to the corresponding OEP values e° EP at similar 

For the Ne atom (Fig. O, the discrep- 



speed as Ap V ' OL 



DV 



DV 



for Ne (Fig. that both indicators A^ 
decrease rapidly with I (like i^^ where C ~ constant) and 
they are of the same order at each step I of the den- 
sity iteration up to I = t\ = 20 where A^y OEP starts 
to saturate at a very small value of the order of 10~ 7 . 
Therefore, a further decrease of A^y'^ -1 ^ for I > l\, 
seen in Fig. 2J and indicating a seemingly better conver- 
gence of nW, is not accompanied by further improving 
of the density (lowering of Aq V ' OEP ) for I > l\. Simi- 
lar ^-dependence of the DV indicators is found for other 
closed-shell atoms. Thus, we conclude that the DV indi- 
cator Apy'^ -1 ^ is a suitable tool to assess the quality of 
the consecutive approximations to n OEP , provided 
that I < l\. In practice, the value of i\ is close to the 
crossing point of the A^y ^ e ~^ vs I line with the Ap gF 
vs I line; see discussion below. With this definition, we 
find i x = 21, 20, 17 for Ne, Ar, Zn. 

In the reported test for the Ne atom, the exchange po- 
tential is calculated using only fc = 2 iteration steps for 
each approximate density n = n^ . Despite so small fco, 
the DSF indicator A^g F decays quickly, like (where 
£i — constant) with increasing I during the density iter- 
ation and it approaches a limiting value, less than 10~ 7 , 
after some number of steps close to i\ (but, a few steps 
earlier for other atoms, like Ar and Zn). The partner in- 
dicator Apg F reaches the same limiting level in the same 
range of £. 

The speed of convergence of iterated density does not 
depend on the type of the ES optimization used in cal- 
culation of v x , as it is seen in Fig. |Ha),(b) for Ne and 
in Table [IT] for Zn. It is also found that increasing the 
number fco of steps in the iteration of the exchange poten- 
tial — although leading to smaller Ap2 F for each density 



W.f- 1 }, A^ } /° EP ancies I 6 "' " e °i EP l decrease with I roughly as ft (i.e. 

with similar exponential decay as A^y' OEP ) and saturate 
at around the 10~ 6 — 10~ 7 Hartree level after a number of 
steps close to l\ = 21. Note that e n i may oscillate around 
e„ EP during iteration, which leads to changing the sign 



of e,, 



,OEP 



- does not accelerates the overall convergence in the 



density iteration. In the test for Zn (Table [TTJ) , the DV 
indicator A^y OEP obtained for fco > 3 is even slightly 
larger, than for fco = 2; it is true for each £ > and 
both types of the ES optimization. Thus, fco = 2 and 
(JV — l)-parameter ES modification is recommended as 
an optimal choice. 

Summarizing the discussion on the use of the DSF and 
DV indicators we conclude that the condition (|6.26[) . with 
suitably chosen accuracy thresholds t/dsf and ?7dVj can 
indeed be used for termination of the combined iteration 



of ra W and 



,[fc] 



In particular, in the discussed test for 



like for (nl) = 2p orbital in the Ne atom. 
The discrepancy of e„; from e^ EP becomes less than 
10 -4 Hartree for all occupied orbitals nl after 9 density 
iterations for Ne, and 12 steps for Ar and Zn while using 
only fco = 2 steps in iteration of v x for each density. Thus, 
the 0.1 mHartree accuracy of orbital energies is obtained 
after a similar number of iterations as in Ref. |28| . where 
a different method of v x iteration is applied and the OS 
are found by solving the KP differential equations. 

The accuracy of calculated energies e n i depends on the 
choice of p and k in the FD formulas (|4.6p . (|4.7[) used 
to determine the approximate ES and OS applied in the 
iteration of t> x . Indeed, in the test for Zn (Table IIII[) . 
the discrepancy of e„; from e° EP decreases rapidly with 
k, especially for p = 3, 4. The discrepancy magnitudes 
depend on k roughly as kP so that they obey the same 
power law as the ES and OS errors corresponding to the 
type of the FD formula applied in the calculations. The 
accuracy of the occupied orbital energies obtained for Zn 
is 10~ 5 and 10~ 6 Hartree for p = 2, k = 0.05 and p = 4, 
k = 0.1, while it reaches 10 -7 Hartree forp = 4, k = 0.05. 
Note that the choice of p, n, though determining the final 
accuracy of e n i , does not affect the rate of convergence of 
the orbital energies during density iteration. 



VIII. CONCLUSIONS 

The performed calculations for closed-(sub)shell atoms 
prove that the exact exchange potential can be deter- 
mined very accurately with the presently proposed algo- 
rithm using only the occupied solutions of the perturbed 
KS equations. Since the respective perturbation is given 
by the (scaled) DFT exchange potential or the (scaled) 
difference between this potential and the non-local Fock 
exchange operator (built with the KS orbitals), the per- 
turbed KS equations have a similar form as the KS or 
HF equations. Thus, these KS-like and HF-like equations 
can be solved with (only slightly modified) standard nu- 
merical methods applicable to the KS and HF equations. 
Therefore, the present algorithm for determination of the 
exact exchange potential, shown in this work to be effec- 
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tive and very accurate in the spatial-grid representation 
(with the pseudospectral method), is plausibly suitable 
for its implementation within standard codes using an 
orbital basis to represent the KS orbitals as well as the 
KS-likc and HF-like orbitals applied in the iteration of 

The applied method of iterating the exchange poten- 
tial, which is the second new element of our algorithm, 
is found to be convergent only if the modification of its 
OS term, performed at even iteration steps, is accompa- 
nied by the optimization of the ES term at odd steps. 
With both ES and OS modifications included, the rapid 
convergence of the obtained DS norm, decreasing below 
10~ 6 in less than 20 iteration steps, proves high efficiency 
of the algorithm and very high accuracy of w x at con- 
vergence (as the DS norm vanishes for exact v x due to 
the OEP equation). Since the iteration of the exchange 
potential can be reduced to only 2 steps (for each den- 
sity) in the self-consistent density calculations within the 
exchange-only KS scheme, without affecting the quick 
convergence of the iterated density, the orbital energies 
of 10~ 6 — 10 _7 Hartree accuracy are obtained for atoms 
with only the total number of 20 x 2 = 40 iterations of 
u x (starting with the KLI approximation). 

Finally, we expect that the proposed algorithm for de- 
termination of the exact exchange potential, successfully 
tested in calculations for atoms in a spatial-grid repre- 
sentation, should also be equally efficient for molecules, 
and that it could be applied within an orbital-basis rep- 
resentation. 
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Appendix A: Reduction of A-dimensional algebraic 
problem 



to (N — l)-dimensional one 

The Af-dimensional system of algebraic equations for 
unknown d , represented by a matrix equation 

A' d' = b' (Al) 

where A is N x N square symmetric matrix {A ac }^2a, 
while d and b are A^ x 1 column matrices {dc}^) 1 and 
j&c}^) 1 , can be rewritten equivalently as a system of the 
following two equations 

A>o do + Aj d = b , (A2) 
A d + A d Q = b , (A3) 

where A is (N — 1) x (N — 1) square matrix {^c}ac==i> 
d, b and Ao are (N — 1) x 1 column matrices {dc}^^ 1 , 
{b c }c=ii {^oolcrJi 1 , while A ao , d a and b are separate 
elements. The solution of the system of Eqs. (IA2[) and 
(|A3|) is found to be 

_ bo-A^A^b 
d ° Aoo-AJA-Ao ' (M) 



d = A -1 (b - A do) • (A5) 

As we see, for its evaluation it is sufficient to determine 
the (N - 1) X (N - 1) matrix A" 1 instead of the N x N 
matrix (A ) 1 , necessary in the case of direct solution of 
the original Eq. (jAll) . And what is more, in our applica- 
tion the matrix A -1 is common for all steps of iterations 
at fixed density, while the matrix (A ) is specific to 
each step separately. 
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\. both 
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(given in the last row). 
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TABLE II. The indicator A^'° EP = ||r^>-n OEP || of the discrepancy between the approximate electron density at the 
£-th iteration step and the exchange-only OEP density n OEP for the Zn atom. The shown results are obtained for each density 
■nS 1 ^ 1 with fco steps in the iteration of the exchange potential and either (N — 1)- or TV-parameter ES optimization. 







A^' OEP (xlO- 3 ) 








(TV — l)-parameter 




TV-parameter 




£ fc = 2 


fco = 3 


fco = 7 fc = 2 


fco = 3 


fco = 7 






3.07390 


3.07390 


3.07390 


3.07390 


3.07390 


3.07390 


1 


1.27749 


1.71563 


1.82028 


1.27749 


1.84049 


1.83820 


2 


0.95047 


1.06142 


1.07960 


0.79210 


1.09074 


1.08965 


3 


0.49044 


0.61554 


0.64348 


0.46845 


0.65009 


0.64967 


5 


0.17837 


0.21939 


0.22823 


0.16534 


0.23049 


0.23052 


10 


0.01355 


0.01659 


0.01724 


0.01242 


0.01739 


0.01742 


20 


0.00049 


0.00060 


0.00060 


0.00049 


0.00060 


0.00060 


30 


0.00047 


0.00057 


0.00057 


0.00047 


0.00057 


0.00057 



TABLE III. The orbital energies t n \ vs. p , k for the Zn atom obtained — after 25 iterations of — with the algorithm of 
Sec. IVII using fco = 2 iteration steps and (TV — l)-parameter ES optimization in determination of the exchange potential 
for each n {e} . The applied OS and ES are found with the FD formulas (|4.6[) . (|4.7I) using the chosen p, k. The shown values 
are discrepancies of e n ; from the orbital energies e° EP (given in the last row) calculated with the exact exchange potential 
w E [n OEP ] determined with the non-iterative method of Ref. [301 . 











e«i ~ e° EP (10- 4 Hartree) 






V 


K 


(nl) = Is 


2s 


3s 4s 2p 


?yp 


M 



1 


0.2 


-10.640 


-9.992 


-12.802 


3.106 


-8.231 


-12.603 


-10.352 


1 


0.1 


-5.192 


-4.714 


-6.001 


1.474 


-3.813 


-5.889 


-4.808 


1 


0.05 


-2.575 


-2.296 


-2.912 


0.719 


-1.842 


-2.853 


-2.323 


2 


0.2 


-0.056 


-1.075 


-1.885 


0.308 


-0.981 


-1.957 


-1.816 


2 


0.1 


-0.019 


-0.268 


-0.466 


0.077 


-0.246 


-0.484 


-0.449 


2 


0.05 


-0.005 


-0.066 


-0.116 


0.019 


-0.061 


-0.121 


-0.112 


3 


0.2 


0.126 


0.661 


0.927 


-0.131 


0.640 


0.963 


0.914 


3 


0.1 


0.010 


0.065 


0.091 


-0.013 


0.062 


0.095 


0.090 


3 


0.05 


0.001 


0.009 


0.011 


-0.002 


0.007 


0.011 


0.010 


4 


0.2 


0.041 


0.170 


0.225 


-0.031 


0.165 


0.233 


0.222 


4 


0.1 


0.002 


0.011 


0.013 


-0.002 


0.010 


0.013 


0.013 


4 


0.05 


-0.0005 


0.0015 


0.0007 


-0.0004 


0.0006 


0.0005 


0.0007 



e%r -345.755720523 -41.714189169 -4.796168733 -0.292805644 -36.742098912 -3.210661901 -0.537803838 
(Hartree) 
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FIG. 1. Approximate radial OS 5J? n ; = Sxni/r = ' p,fv 'xni/ r calculated for the Ne atom with the FD formulas for various p, n 
and compared with the exact OEP OS 5i?° EP , all obtained for the the exchange-only OEP density n OEP and the corresponding 
exact exchange potential u E [n OEP ] determined with the non-iterative method of Ref. [30L 
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FIG. 2. (a,b) Density shift <5n' fe ' for the approximate exchange potential i4 fc ' , (c) difference between t4 fe ' and the exact exchange 
potential , all obtained for the the exchange-only OEP density n OEP in the Ne atom. The plots are shown for k — 0, 2, 4, 6, 8 
(solid lines) and fc = 1,3,5,7 (dotted lines) and each of them is labeled with the corresponding k value. For even k, the 
respective is obtained by modifying the OS term in li* -1 ' while for odd k it is found by (N — l)-parameter optimization 
of the ES of v l J" 1] . For better comparison with an additional constant shift is applied to «I fe ' at each iteration step k to 



obtain D 



0, the same as Di\ 
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10 20 

k (iterations of v x ) 

FIG. 3. DSF indicator A^g F = IcW^H vs. the index k of the iterative approximation ti*' to the exact exchange potential for 
the exchange-only OEP density n OEP in (a) Ne, (b) Ar and (c) Zn atoms. The respective is determined with the algorithm 
of Sec. IVII for even k — by modifying the OS term in for odd k — by the ES term modification: the (N — l)-parameter 

optimization of the ES (solid circles), the JV-parameter optimization (open circles), no optimization, i.e., empty step (open 
triangles). 
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FIG. 4. DSF indicators A l ° ] SF = ||<5n [01 || (open circles), = \\Sn lk ° ] \\ (solid circles) and DV indicators A 



{f },OEP 
DV 



i HEP || (solid squares), A^'" 1 ^ = \\n^ — n^ 1 ^ | (open squares) for consecutive approximations to the exchange- 



only OEP electron density n OEP in the Ne atom. The results are obtained for the approximate exchange potentials found — 
for each — with the algorithm of Sec. IVII using ko = 2 steps in the iteration of 4*' and employing (a) the (N — l)-parameter 
optimization of the ES, and (b) the iV-parameter optimization. The starting density n<°> is found with the self-consistent KLI 
exchange potential • The dotted line follows the path of the the u x iteration. 
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FIG. 5. The absolute differences between the approximate orbital energies e n i obtained for the Ne atom with 
the KS potential vi^ , Eq. (|6.25[) . at the ^-th density iteration step, and the respective energies e° EP (equal to 
-30.8200039329, -1.7181258001, -0.8507101622 Hartree for (nl) = ls,2s,2p) found with the exact exchange potential for 
n OEP . The exchange potential 4*°' entering is found (for n^ -1 ^) with ko — 2 steps of iteration which employs the 
(N — l)-parameter optimization of the ES. The results for nl — Is (circles), 2s (squares) and 2p (triangles) are marked with 
solid symbols for positive values of e„j — eS EP and with open symbols for negative ones. 



